LPA and SC

suppressPackageStartupMessages({
  library("here")
  library("tidyverse")
  library("MplusAutomation")
  library("gt")
  library("glue")
  library("kableExtra")
  library("misty")
  library("lavaan")
  library("AICcmodavg")
  library("nonnest2")
  library("DiagrammeR")
  library("lavaan")
  library("tidyLPA")
  library("semTools")
  library("brms")
  library("MBESS")
  library("ufs")
  library("robmed")
  library("careless")
  library("psych")
  library("BayesFactor")
  library("effectsize")
  library("tidybayes")
  library("emmeans")
  library("bayesplot")
  library("patchwork")
  library("bmlm")
})

options("max.print" = .Machine$integer.max)

# Make random things reproducible
set.seed(1234)

options(mc.cores = 4)

bayesplot_theme_set()

source(here::here("src", "R", "functions", "funs_add_neoffi60_subscales.R"))
source(here::here("src", "R", "functions", "funs_correct_iesr_scores.R"))
source(here::here("src", "R", "functions", "funs_plot_job_qualification.R"))
source(here::here("src", "R", "functions", "funs_generate_all_items_df.R"))

scale_this <- function(x) as.vector(scale(x))

sum_coding <- function(x, lvls = levels(x)) {
  # codes the first category with -1
  nlvls <- length(lvls)
  stopifnot(nlvls > 1)
  cont <- diag(nlvls)[, -nlvls, drop = FALSE]
  cont[nlvls, ] <- -1
  cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
  colnames(cont) <- lvls[-1]
  x <- factor(x, levels = lvls)
  contrasts(x) <- cont
  x
}

Get data

all_items <- generate_all_items_df()

There is a problem with IES-R, in the control group. I shift the control distribution of IES-R towards lower values.

temp <- correct_iesr_scores(all_items)
all_items <- temp

ggplot(all_items, aes(x = iesr_ts, colour = is_rescue_worker)) +
  geom_density()

all_items |> 
  group_by(is_rescue_worker) |> 
  summarize(
    avg_iesr = mean(iesr_ts)
  )

There are some cases in which the RWs of Toscana and Lombardia did not have the proper qualification. They are assigned to the category of RWs.

all_items$idx <- 1:nrow(all_items)
all_items$is_rescue_worker <- as.character(all_items$is_rescue_worker)
all_items$is_rescue_worker <- ifelse(
  all_items$idx < 746, "Si", all_items$is_rescue_worker
)
all_items$idx <- NULL

all_items$commeetee_location <- ifelse(
  all_items$is_rescue_worker == "No", "None", all_items$red_cross_commeetee_location
) |> 
  factor()

all_items$commeetee_location <- ifelse(
  all_items$is_rescue_worker == "No" & all_items$age < 25, "students", all_items$commeetee_location
)

all_items$commeetee <- ifelse(
  all_items$is_rescue_worker == "No" & all_items$age >= 25, "community_sample",
  all_items$commeetee_location
)

all_items$commeetee <- 
  ifelse(is.na(all_items$commeetee), "community_sample", all_items$commeetee) |> 
  factor()

sum(is.na(all_items$commeetee))

[1] 0

# all_items <- all_items %>%
#   mutate(job_qualification = case_when(
#     is_rescue_worker=="Si" & job_qualification == "non_rescue_worker" ~ "team_member",
#     TRUE ~ job_qualification)) 

IES-R as a function of group

iesr_ts | trunc(lb = 0) ~ is_rescue_worker + (1 | commeetee),

m0 <- brm(
  bf(
    iesr_ts ~ is_rescue_worker,
    sigma ~ is_rescue_worker 
    ),
  family = skew_normal(),
  data = all_items,
  backend = "cmdstanr"
)

Compiling Stan program...
Warning: incomplete final line found on '/var/folders/hl/dt523djx7_q7xjrthzjpdvc40000gn/T//Rtmp9PpTYj/model-122363e511137.hpp'

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

# m0 <- brm(
#   bf(
#     iesr_ts ~ is_rescue_worker + (1 | commeetee),
#     sigma ~ is_rescue_worker + (1 | commeetee)
#     ),
#   family = skew_normal(),
#   data = all_items,
#   backend = "cmdstanr"
# )

pp_check(m0)

Error: object 'm0' not found

summary(m0)

 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: iesr_ts ~ is_rescue_worker 
         sigma ~ is_rescue_worker
   Data: all_items (Number of observations: 1068) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                   11.72      0.50    10.80    12.78 1.00     1941     2352
sigma_Intercept              2.31      0.04     2.23     2.39 1.00     1966     2109
is_rescue_workerSi           7.16      0.70     5.80     8.54 1.00     1706     1949
sigma_is_rescue_workerSi     0.37      0.05     0.27     0.46 1.00     1828     2083

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
alpha    21.73      2.40    17.27    26.63 1.00     2632     2121

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m0, "is_rescue_worker"
)
plot(me, points = FALSE)

BFt <- BayesFactor::ttestBF(
  all_items$ies_ts[all_items$is_rescue_worker == "Si"], 
  all_items$ies_ts[all_items$is_rescue_worker == "No"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.69      | [0.55, 0.82]

- Estimated using pooled SD.

Supported families are: ‘acat’, ‘asym_laplace’, ‘bernoulli’, ‘beta’, ‘beta_binomial’, ‘binomial’, ‘categorical’, ‘com_poisson’, ‘cox’, ‘cratio’, ‘cumulative’, ‘custom’, ‘dirichlet’, ‘dirichlet2’, ‘discrete_weibull’, ‘exgaussian’, ‘exponential’, ‘frechet’, ‘gamma’, ‘gaussian’, ‘gen_extreme_value’, ‘geometric’, ‘hurdle_cumulative’, ‘hurdle_gamma’, ‘hurdle_lognormal’, ‘hurdle_negbinomial’, ‘hurdle_poisson’, ‘info’, ‘inverse.gaussian’, ‘logistic_normal’, ‘lognormal’, ‘multinomial’, ‘negbinomial’, ‘negbinomial2’, ‘poisson’, ‘shifted_lognormal’, ‘skew_normal’, ‘sratio’, ‘student’, ‘von_mises’, ‘weibull’, ‘wiener’, ‘zero_inflated_asym_laplace’, ‘zero_inflated_beta’, ‘zero_inflated_beta_binomial’, ‘zero_inflated_binomial’, ‘zero_inflated_negbinomial’, ‘zero_inflated_poisson’, ‘zero_one_inflated_beta’

The sk, ch, mi sub-scales are coded so that high values indicate high self-compassion levels. The sj, is, oi sub-scales are coded so that high values indicate low self-compassion levels.

The ts_sc score has been computed by reversing the coding of the items of the sj, is, oi sub-scales (so that they indicate the absence of self-judgment, absence of isolation, absence of over-identification).

scs_subscales <- with(all_items, data.frame(sk, ch, mi, sj, is, oi, scs_ts))
cor(scs_subscales) |> round(2)

          sk    ch    mi    sj    is    oi scs_ts
sk      1.00  0.52  0.58 -0.39 -0.28 -0.24   0.71
ch      0.52  1.00  0.49 -0.01 -0.03 -0.04   0.45
mi      0.58  0.49  1.00 -0.19 -0.33 -0.35   0.66
sj     -0.39 -0.01 -0.19  1.00  0.67  0.66  -0.75
is     -0.28 -0.03 -0.33  0.67  1.00  0.80  -0.78
oi     -0.24 -0.04 -0.35  0.66  0.80  1.00  -0.78
scs_ts  0.71  0.45  0.66 -0.75 -0.78 -0.78   1.00

COPE scale

In the COPE scale only two factors are identified.

all_items$pos_reinterpretation <- with(all_items, cope_1 + cope_29 + cope_38 + cope_59)
all_items$mental_disengagement <- with(all_items, cope_2 + cope_16 + cope_31 + cope_43) 
all_items$venting <- with(all_items, cope_3 + cope_17 + cope_28 + cope_46) 
all_items$seeking_instrumental_support <- with(all_items, cope_4 + cope_14 + cope_30 + cope_45) 
all_items$active_coping <- with(all_items, cope_5 + cope_25 + cope_47 + cope_58)  
all_items$denial <- with(all_items, cope_6 + cope_27 + cope_40 + cope_57) 
all_items$religion <- with(all_items, cope_7 + cope_18 + cope_48 + cope_60) 
all_items$humor <- with(all_items, cope_8 + cope_20 + cope_36 + cope_50) 
all_items$behavioral_disengagement <- with(all_items, cope_9 + cope_24 + cope_37 + cope_51) 
all_items$restraint <- with(all_items, cope_10 + cope_22 + cope_41 + cope_49) 
all_items$seeking_emotional_support <- with(all_items, cope_11 + cope_23 + cope_34 + cope_52) 
all_items$substance_use <- with(all_items, cope_12 + cope_26 + cope_35 + cope_53) 
all_items$acceptance <- with(all_items, cope_13 + cope_21 + cope_44 + cope_54) 
all_items$suppr_competing_activities <- with(all_items, cope_15 + cope_33 + cope_42 + cope_55) 
all_items$planning <- with(all_items, cope_19 + cope_32 + cope_39 + cope_56) 

Create COPE sub-scales scores using all items – note that SEM analyses suggest to drop some of the items.

all_items$active_coping <- with(
  all_items, pos_reinterpretation + active_coping +
  suppr_competing_activities + planning + restraint + 
    seeking_instrumental_support + acceptance
)

all_items$avoidance_coping <- with(
  all_items, mental_disengagement + denial + humor +
  behavioral_disengagement + substance_use + religion 
)

all_items$soc_emo_coping <- with(
  all_items, seeking_instrumental_support +
  seeking_emotional_support + venting
)

Self-compassion scale

plot(density(all_items$scs_ts))

pp_check(fit_1)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

me <- conditional_effects(
  fit_1, "is_rescue_worker"
)
plot(me, points = FALSE)

summary(fit_1)

 Family: student 
  Links: mu = identity; sigma = log; nu = identity 
Formula: scs_ts ~ is_rescue_worker 
         sigma ~ is_rescue_worker
   Data: all_items (Number of observations: 1068) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                   74.13      1.01    72.17    76.08 1.00     4836
sigma_Intercept              2.88      0.04     2.81     2.97 1.00     4420
is_rescue_workerSi           7.45      1.16     5.15     9.67 1.00     4717
sigma_is_rescue_workerSi    -0.13      0.05    -0.22    -0.03 1.00     5247
                         Tail_ESS
Intercept                    2949
sigma_Intercept              2748
is_rescue_workerSi           2725
sigma_is_rescue_workerSi     2970

Family Specific Parameters: 
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu    38.43     15.26    16.30    75.45 1.00     4067     3086

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

BFt <- BayesFactor::ttestBF(
  all_items$scs_ts[all_items$is_rescue_worker == "Si"], 
  all_items$scs_ts[all_items$is_rescue_worker == "No"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.42      | [0.29, 0.56]

- Estimated using pooled SD.

rw_df <- all_items |> 
  dplyr::filter(is_rescue_worker == "Si")

rw_df <- rw_df %>%
  mutate(job_qualification = case_when(
    job_qualification == "non_rescue_worker" ~ "team_member",
    TRUE ~ job_qualification))

LPA

lpa_scales <- c(
  "is_rescue_worker",
  "neuroticism", "extraversion", "openness", "agreeableness", "conscientiousness",
  "active_coping", "avoidance_coping", "soc_emo_coping",
  "iesr_ts",
  # "avoiding", "intrusivity", "hyperarousal",
  # "sk", "ch", "mi", "sj", "is", "oi",
  # "pos_sc",
  # "neg_sc",
  # "ts_sc",
  "mpss_tot"
  # "ptgi_total_score"
  # "relating_to_others",
  # "new_possibilities",
  # "personal_strength",
  # "appreciation_of_life",
  # "spirituality"
)

lpa_rw_df <- subset(rw_df, select=lpa_scales) 

# lpa_rw_df <- subset(all_items, select=lpa_scales) |> 
#   dplyr::filter(is_rescue_worker == "Si")

lpa_rw_df$is_rescue_worker <- NULL

```r
lpa_rw_df %>% 
  scale() %>%
  estimate_profiles(1:6,
    variances = c(\equal\, \varying\),
    covariances = c(\zero\, \varying\)
    #package = \MplusAutomation\
  )

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


 Model Classes AIC      BIC      Entropy prob_min prob_max n_min n_max BLRT_p
 1     1       21200.56 21292.85 1.00    1.00     1.00     1.00  1.00        
 1     2       20604.04 20747.10 0.71    0.87     0.94     0.36  0.64  0.00  
 1     3       20322.46 20516.27 0.75    0.82     0.93     0.18  0.56  0.00  
 1     4       20180.65 20425.23 0.74    0.83     0.86     0.10  0.42  0.00  
 1     5       20066.42 20361.76 0.74    0.79     0.88     0.08  0.34  0.00  
 1     6       20007.30 20353.40 0.73    0.76     0.86     0.06  0.30  0.00  
 1     7       19941.06 20337.92 0.76    0.77     0.85     0.03  0.31  0.00  
 1     8       19877.85 20325.48 0.78    0.74     0.88     0.02  0.30  0.00  
 1     9       19826.21 20324.60 0.78    0.75     0.91     0.02  0.33  0.00  
 1     10      19793.42 20342.58 0.78    0.75     0.92     0.01  0.33  0.00  
 6     1       19747.98 20047.94 1.00    1.00     1.00     1.00  1.00        
 6     2       19380.08 19984.61 0.70    0.89     0.94     0.45  0.55  0.00  
 6     3       19286.96 20196.06 0.74    0.88     0.89     0.23  0.39  0.00  
 6     4       19262.78 20476.45 0.78    0.83     0.93     0.11  0.43  1.00  
 6     5       19282.17 20800.41 0.78    0.81     0.92     0.07  0.43  0.43  
 6     6       19257.62 21080.44 0.86    0.83     1.00     0.03  0.51  0.00  
 6     7       19211.10 21338.49 0.84    0.78     1.00     0.01  0.47        
 6     8       19629.33 22061.29 0.86    0.00     0.95     0.00  0.53        
 6     9                                                                     
 6     10      19648.29 22675.55 0.87    0.00     1.00     0.00  0.42     


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->



<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Compare tidyLPA solutions:

 Model Classes AIC       BIC      
 1     1       21342.450 21434.878
 1     2       20746.522 20889.786
 1     3       20464.673 20658.772
 1     4       20341.265 20586.199
 1     5       20245.437 20541.207
 1     6       20175.056 20521.661
 6     1       19889.997 20190.388
 6     2       19515.286 20120.691
 6     3       19437.181 20347.598
 6     4       19389.596 20605.025
 6     5       19400.107 20920.550
 6     6       19404.469 21229.925

Best model according to AIC is Model 6 with 4 classes.
Best model according to BIC is Model 6 with 2 classes.

An analytic hierarchy process, based on the fit indices AIC, AWE, BIC, CLC, and KIC (Akogul & Erisoglu, 2017), suggests the best solution is Model 6 with 2 classes.



<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubTIgPC0gbHBhX3J3X2RmICU+JVxuICBzY2FsZSgpICU+JVxuICBlc3RpbWF0ZV9wcm9maWxlcygyLFxuICAgIHZhcmlhbmNlcyA9IFwidmFyeWluZ1wiLFxuICAgIGNvdmFyaWFuY2VzID0gXCJ2YXJ5aW5nXCIsXG4gICAgcGFja2FnZSA9IFwiTXBsdXNBdXRvbWF0aW9uXCJcbiAgKVxuYGBgIn0= -->

```r
m2 <- lpa_rw_df %>%
  scale() %>%
  estimate_profiles(2,
    variances = "varying",
    covariances = "varying",
    package = "MplusAutomation"
  )

m2_plot <- lpa_rw_df %>%
  scale() %>%
  estimate_profiles(2,
    variances = "varying",
    covariances = "varying"
    #package = "MplusAutomation"
    ) %>%
    plot_profiles(add_line = TRUE, rawdata= FALSE, bw = FALSE)

Profile 1: dysfunctional Profile 2: adaptive

get_estimates(m2)

out <- get_data(m2)
lpa_rw_df$lpa_class <- out$Class

table(
  lpa_rw_df$lpa_class
)

  1   2 
332 419 

table(
  lpa_rw_df$lpa_class, rw_df$job_qualification
)

   
    driver team_leader team_member
  1     57         135         140
  2     98         164         157

lpa_rw_df$class <- factor(lpa_rw_df$lpa_class)
summary(lpa_rw_df$class)

  1   2 
332 419 

rw_df$class <- lpa_rw_df$class

m1 <- brm(
  bf(scs_ts ~ class),
  family = student(),
  data = rw_df,
  init = 0.1,
  backend = "cmdstanr",
  adapt_delta = 0.9
)

pp_check(m1)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m1)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: scs_ts ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    74.19      0.83    72.60    75.82 1.00     3812     3100
class2       13.07      1.10    11.00    15.21 1.00     4001     3053

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    14.30      0.47    13.38    15.20 1.00     3617     2817
nu       28.45     13.18    11.35    62.77 1.00     3428     2922

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m1, "class"
)
plot(me, points = FALSE)

BFt <- BayesFactor::ttestBF(
  rw_df$scs_ts[rw_df$class == 1], 
  rw_df$scs_ts[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt, type = "d")

Cohen's d |         95% CI
--------------------------
-0.87     | [-1.02, -0.72]

- Estimated using pooled SD.

m1 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  labs(x = "LPA Class", y = "SCS Score", title = "Rescue Workers") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 78.5, label = "Bayesian Cohen's d = 0.89\n 95% CI [0.73, 1.04]")

```r
names(all_items)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


### Self-judgment


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubTIgPC0gYnJtKFxuICBiZihzaiB+IGNsYXNzKSxcbiAgZGF0YSA9IHJ3X2RmLFxuICBmYW1pbHkgPSBzdHVkZW50LFxuICBiYWNrZW5kID0gXCJjbWRzdGFuclwiLFxuICBhZGFwdF9kZWx0YSA9IDAuOTlcbilcbmBgYCJ9 -->

```r
m2 <- brm(
  bf(sj ~ class),
  data = rw_df,
  family = student,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m2)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m2)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: sj ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    17.12      0.25    16.63    17.62 1.00     3494     2516
class2       -3.67      0.32    -4.29    -3.06 1.00     3404     2787

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.26      0.12     4.03     4.51 1.00     3787     2427
nu       46.52     16.90    21.85    86.76 1.00     3495     2719

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m2, "class"
)
plot(me, points = FALSE)

emmeans(m2, specs = pairwise ~ class)

$emmeans
 class emmean lower.HPD upper.HPD
 1       17.1      16.6      17.6
 2       13.4      13.0      13.9

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast        estimate lower.HPD upper.HPD
 class1 - class2     3.66      3.07       4.3

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$sj[rw_df$class == 1], 
  rw_df$sj[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.83      | [0.67, 0.98]

- Estimated using pooled SD.

p2 <- m2 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Self-Judgment", title = "A") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 16, label = "Bayesian Cohen's d = 0.82\n 95% CI [0.67, 0.97]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p2

Isolation

m3 <- brm(
  bf(is ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m3)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m3)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: is ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    12.37      0.25    11.86    12.87 1.00     2367     2295
class2       -3.45      0.33    -4.11    -2.80 1.00     2399     2641

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.05      0.12     3.83     4.30 1.00     2133     2366
alpha     2.19      0.66     0.88     3.46 1.00     1966     1683

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m3, "class"
)
plot(me, points = FALSE)

emmeans(m3, specs = pairwise ~ class)

$emmeans
 class emmean lower.HPD upper.HPD
 1      12.37     11.84      12.8
 2       8.93      8.54       9.3

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast        estimate lower.HPD upper.HPD
 class1 - class2     3.44      2.81      4.11

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$is[rw_df$class == 1], 
  rw_df$is[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.94      | [0.79, 1.10]

- Estimated using pooled SD.

p3 <- m3 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Isolation", title = "B") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 11, label = "Bayesian Cohen's d = 0.94\n 95% CI [0.79, 1.10]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p3

Over-identification

m4 <- brm(
  bf(oi ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m4)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m4)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: oi ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    11.10      0.25    10.58    11.57 1.00     1434     1663
class2       -2.68      0.35    -3.35    -1.95 1.01     1372     1365

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.63      0.12     3.41     3.87 1.00     1445     1791
alpha     3.24      0.91     1.94     5.48 1.01     1314     1122

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m4, "class"
)
plot(me, points = FALSE)

emmeans(m4, specs = pairwise ~ class)

$emmeans
 class emmean lower.HPD upper.HPD
 1      11.11     10.60     11.58
 2       8.42      8.05      8.79

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast        estimate lower.HPD upper.HPD
 class1 - class2     2.69      2.01       3.4

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$oi[rw_df$class == 1], 
  rw_df$oi[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.99      | [0.83, 1.15]

- Estimated using pooled SD.

p4 <- m4 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('High Resilience', 'Low Resilience')) +
  labs(x = "LPA Class", y = "SCS Over-Identification", title = "C") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 9.4, label = "Bayesian Cohen's d = 0.97\n 95% CI [0.82, 1.12]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p4

Self-kindness

m5 <- brm(
  bf(sk ~ class),
  family = student(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m5)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m5)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: sk ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    13.07      0.24    12.59    13.54 1.00     3482     2388
class2        1.18      0.32     0.55     1.80 1.00     3684     2519

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.22      0.11     4.00     4.45 1.00     3537     2868
nu       41.60     16.59    17.38    80.94 1.00     3185     2731

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m5, "class"
)
plot(me, points = FALSE)

emmeans(m5, specs = pairwise ~ class)

$emmeans
 class emmean lower.HPD upper.HPD
 1       13.1      12.6      13.5
 2       14.2      13.9      14.7

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast        estimate lower.HPD upper.HPD
 class1 - class2    -1.17     -1.79    -0.545

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$sk[rw_df$class == 1], 
  rw_df$sk[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |         95% CI
--------------------------
-0.27     | [-0.41, -0.13]

- Estimated using pooled SD.

p5 <- m5 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('High Resilience', 'Low Resilience')) +
  labs(x = "LPA Class", y = "SCS Self-Kindness", title = "D") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 14.5, label = "Bayesian Cohen's d = 0.28\n 95% CI [0.14, 0.43]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p5

Common humanity

m6 <- brm(
  bf(ch ~ class),
  family = student(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m6)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m6)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: ch ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    11.61      0.18    11.27    11.96 1.00     3502     2643
class2       -0.04      0.25    -0.53     0.45 1.00     3016     2397

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.28      0.09     3.10     3.47 1.00     2680     2104
nu       42.34     16.50    18.20    80.94 1.00     2954     2163

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m6, "class"
)
plot(me, points = FALSE)

BFt <- BayesFactor::ttestBF(
  rw_df$ch[rw_df$class == 1], 
  rw_df$ch[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |        95% CI
-------------------------
9.51e-03  | [-0.13, 0.15]

- Estimated using pooled SD.

p6 <- m6 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('High Resilience', 'Low Resilience')) +
  labs(x = "LPA Class", y = "SCS Common-Humanity", title = "E") +
  # papaja::theme_apa() + 
  annotate("text", x = 1.5, y = 12.3, label = "Bayesian Cohen's d = 0.00\n 95% CI [-0.14, 0.14]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p6

Mindfulness

m7 <- brm(
  bf(mi ~ class),
  family = student(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m7)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m7)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: mi ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    12.76      0.17    12.43    13.09 1.00     3277     2803
class2        0.96      0.23     0.51     1.42 1.00     3209     2387

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.04      0.09     2.87     3.22 1.00     2578     2822
nu       36.60     15.48    15.14    74.03 1.00     2806     3032

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

emmeans(m7, specs = pairwise ~ class)

$emmeans
 class emmean lower.HPD upper.HPD
 1       12.8      12.4      13.1
 2       13.7      13.4      14.0

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast        estimate lower.HPD upper.HPD
 class1 - class2   -0.964     -1.43     -0.52

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$mi[rw_df$class == 1], 
  rw_df$mi[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |         95% CI
--------------------------
-0.30     | [-0.44, -0.15]

- Estimated using pooled SD.

me <- conditional_effects(
  m7, "class"
)
plot(me, points = FALSE)

p7 <- m7 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Mindfulness", title = "F") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 13.8, label = "Bayesian Cohen's d = 0.32\n 95% CI [0.17, 0.46]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p7

fig_scs <- (p2 | p3 | p4) /
(p5 | p6 | p7)

out <- fig_scs + plot_annotation(
  title = 'SCS Subscales as a Function of LPA Class',
  subtitle = 'Rescue Workers group'
  # caption = 'Disclaimer: None of these plots are insightful'
)
ggsave("scs_subscales_lpa.pdf", width = 35, height = 20, units = "cm")

IES-R

m10 <- brm(
  bf(ies_ts ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m10)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m10)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: ies_ts ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    21.08      0.72    19.77    22.68 1.00     1705     1915
class2       -4.26      1.05    -6.71    -2.57 1.00     1211     1474

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    13.95      0.40    13.19    14.75 1.00     1229     1470
alpha    10.07      2.34     5.90    15.14 1.00      973     1426

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

p10 <- m10 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "Impact of Event Scale - Revised (IES-R)") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 17, label = "Bayesian Cohen's d = 1.34\n 95% CI [1.18, 1.50]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p10

BFt <- BayesFactor::ttestBF(
  rw_df$ies_ts[rw_df$class == 1], 
  rw_df$ies_ts[rw_df$class == 2],
  paired = FALSE
)

t is large; approximation invoked.

effectsize(BFt)

Cohen's d |       95% CI
------------------------
1.37      | [1.22, 1.53]

- Estimated using pooled SD.

emmeans(m10 , specs = pairwise ~ class)

$emmeans
 class emmean lower.HPD upper.HPD
 1       21.0      19.6      22.4
 2       16.8      15.4      18.2

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast        estimate lower.HPD upper.HPD
 class1 - class2     4.13      2.38      6.41

Point estimate displayed: median 
HPD interval probability: 0.95 

```r
m11 <- brm(
  bf(ptgi_total_score | trunc(lb = 0) ~ class),
  family = student(),
  data = rw_df,
  backend = \cmdstanr\
)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHBfY2hlY2sobTExKVxuYGBgXG5gYGAifQ== -->

```r
```r
pp_check(m11)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc3VtbWFyeShtMTEpXG5gYGBcbmBgYCJ9 -->

```r
```r
summary(m11)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


### Job qualification


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucndfZGYkam9iX3F1YWxpZmljYXRpb24gPC0gaWZlbHNlKFxuICByd19kZiRqb2JfcXVhbGlmaWNhdGlvbiA9PSBcIm5vbl9yZXNjdWVfd29ya2VyXCIsIFwidGVhbV9tZW1iZXJcIiwgXG4gIHJ3X2RmJGpvYl9xdWFsaWZpY2F0aW9uXG4pIFxuYGBgIn0= -->

```r
rw_df$job_qualification <- ifelse(
  rw_df$job_qualification == "non_rescue_worker", "team_member", 
  rw_df$job_qualification
) 

m12 <- brm(
  bf(ies_ts ~ job_qualification),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m12)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m12)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: ies_ts ~ job_qualification 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                             Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept                       17.91      0.68    16.56    19.24 1.00
job_qualificationteam_leader     1.20      0.65    -0.02     2.49 1.00
job_qualificationteam_member     1.35      0.66     0.12     2.68 1.00
                             Bulk_ESS Tail_ESS
Intercept                        1914     2193
job_qualificationteam_leader     2133     1952
job_qualificationteam_member     2215     2001

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    14.47      0.39    13.74    15.28 1.00     1942     2043
alpha    15.62      2.23    11.69    20.44 1.00     2406     2098

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m12, "job_qualification"
)
plot(me, points = FALSE)

emmeans(m12 , specs = pairwise ~ job_qualification)

$emmeans
 job_qualification emmean lower.HPD upper.HPD
 driver              17.9      16.5      19.2
 team_leader         19.1      18.0      20.3
 team_member         19.2      18.1      20.4

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast                  estimate lower.HPD upper.HPD
 driver - team_leader        -1.182     -2.45    0.0631
 driver - team_member        -1.345     -2.67   -0.1136
 team_leader - team_member   -0.149     -1.27    0.8995

Point estimate displayed: median 
HPD interval probability: 0.95 

m13 <- brm(
  bf(scs_ts ~ job_qualification * class),
  family = gaussian(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m13)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

me <- conditional_effects(
  m13, "job_qualification:class"
)
plot(me, points = FALSE)

summary(m13)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: scs_ts ~ job_qualification * class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                                    Estimate Est.Error l-95% CI u-95% CI
Intercept                              78.77      1.90    75.10    82.59
job_qualificationteam_leader           -3.78      2.30    -8.45     0.61
job_qualificationteam_member           -7.34      2.26   -11.84    -3.04
class2                                 10.15      2.42     5.40    14.80
job_qualificationteam_leader:class2     2.03      2.99    -3.76     7.96
job_qualificationteam_member:class2     4.48      2.94    -1.42    10.19
                                    Rhat Bulk_ESS Tail_ESS
Intercept                           1.00     1156     1812
job_qualificationteam_leader        1.00     1185     1852
job_qualificationteam_member        1.00     1235     1860
class2                              1.00     1080     1717
job_qualificationteam_leader:class2 1.00     1136     1696
job_qualificationteam_member:class2 1.00     1156     1799

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    14.80      0.39    14.06    15.59 1.00     2650     2594

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

emmeans(m13 , specs = pairwise ~ job_qualification*class)

$emmeans
 job_qualification class emmean lower.HPD upper.HPD
 driver            1       78.7      75.2      82.6
 team_leader       1       75.0      72.4      77.3
 team_member       1       71.4      69.0      73.8
 driver            2       88.9      86.2      91.9
 team_leader       2       87.2      84.9      89.4
 team_member       2       86.1      83.8      88.3

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast                                estimate lower.HPD upper.HPD
 driver class1 - team_leader class1          3.72    -0.306      8.71
 driver class1 - team_member class1          7.27     3.137     11.88
 driver class1 - driver class2             -10.24   -14.752     -5.37
 driver class1 - team_leader class2         -8.43   -12.880     -4.12
 driver class1 - team_member class2         -7.34   -11.794     -3.05
 team_leader class1 - team_member class1     3.58     0.088      7.06
 team_leader class1 - driver class2        -13.93   -18.127    -10.51
 team_leader class1 - team_leader class2   -12.17   -15.285     -8.69
 team_leader class1 - team_member class2   -11.06   -14.490     -7.82
 team_member class1 - driver class2        -17.50   -21.106    -13.43
 team_member class1 - team_leader class2   -15.79   -19.155    -12.51
 team_member class1 - team_member class2   -14.67   -17.974    -11.21
 driver class2 - team_leader class2          1.75    -1.786      5.66
 driver class2 - team_member class2          2.85    -1.065      6.51
 team_leader class2 - team_member class2     1.12    -1.997      4.29

Point estimate displayed: median 
HPD interval probability: 0.95 

m14 <- brm(
  bf(ptgi_total_score ~ ies_total_score | (1 | comme\)),
  family = gaussian(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

summary(m14)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: ptgi_total_score ~ ies_total_score 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          33.36      1.30    30.82    35.89 1.00     2784     2099
ies_total_score     0.40      0.05     0.29     0.50 1.00     2814     2352

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    23.07      0.61    21.89    24.29 1.00     3030     2390

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

bayes_R2(m14)

    Estimate  Est.Error       Q2.5     Q97.5
R2 0.1089372 0.01913069 0.07305561 0.1471529

conditional_effects(m14, "is:class")

SEM

mydf <- data.frame(
  scs = scale(rw_df$scs_ts),
  class = ifelse(rw_df$class == 1, 0.0, 1.0),
  ptgi = scale(rw_df$ptgi_total_score),
  psc = scale(rw_df$sk + rw_df$ch + rw_df$mi),
  nsc = scale(rw_df$sj + rw_df$oi + rw_df$is),
  commettee = rw_df$red_cross_commeetee_location
)

mydf <- mydf[complete.cases(mydf), ]

# Impute NAs with the mode.
# mydf$commettee[is.na(mydf$commettee)] <- "Comitato di Groane"
# summary(factor(mydf$commettee))

f1 <- bf(scs ~ class, family = skew_normal())
f2 <- bf(ptgi ~ scs + class, family = skew_normal())
mod <- brm(
  f1 + f2 + set_rescor(FALSE), 
  data = mydf, 
  cores = 4, 
  refresh = 0,
  init = 0.1,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

bayestestR::mediation(mod, mediator = "scs", ci = 0.95, method = "SPI")

# Causal Mediation Analysis for Stan Model

  Treatment: class
  Mediator : scs
  Response : ptgi

Effect                 | Estimate |          95% SPI
----------------------------------------------------
Direct Effect (ADE)    |   -0.303 | [-0.461, -0.144]
Indirect Effect (ACME) |    0.120 | [ 0.055,  0.184]
Mediator Effect        |    0.153 | [ 0.070,  0.228]
Total Effect           |   -0.183 | [-0.337, -0.049]

Proportion mediated: -65.54% [-257.06%, 125.97%]

Direct and indirect effects have opposite directions. The proportion mediated is not meaningful.

pp_check(mod, resp = "ptgi")

Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check(mod, resp = "scs")

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(mod)

 Family: MV(skew_normal, skew_normal) 
  Links: mu = identity; sigma = identity; alpha = identity
         mu = identity; sigma = identity; alpha = identity 
Formula: scs ~ class 
         ptgi ~ scs + class 
   Data: mydf (Number of observations: 746) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
scs_Intercept     -0.44      0.05    -0.54    -0.35 1.00     3911     2852
ptgi_Intercept     0.17      0.06     0.05     0.29 1.00     3619     3231
scs_class          0.79      0.07     0.65     0.92 1.00     3666     2804
ptgi_scs           0.15      0.04     0.07     0.23 1.00     3199     2879
ptgi_class        -0.30      0.08    -0.47    -0.15 1.00     3176     2768

Family Specific Parameters: 
           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_scs      0.91      0.02     0.87     0.96 1.00     4070     2612
sigma_ptgi     0.99      0.03     0.94     1.04 1.00     4207     2703
alpha_scs     -0.70      0.59    -1.52     0.59 1.00     2244     3333
alpha_ptgi     0.09      0.53    -0.88     1.10 1.00     3801     3079

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

fit <- mlm(d = mydf, 
           id = "commettee",
           x = "class",
           m = "scs",
           y = "ptgi",
           iter = 2000, 
           cores = 4)

mlm_path_plot(fit, level = .95, text = T,
              xlab = "Resilience\nProfile",
              mlab = "Self\nCompassion",
              ylab = "PTG", digits = 2)

Causal Mediation Analysis for Stan Model

Treatment: class Mediator : scs Response : ptgi

Effect | Estimate | 95% SPI

Direct Effect (ADE) | -0.300 | [-0.464, -0.140] Indirect Effect (ACME) | 0.121 | [ 0.053, 0.187] Mediator Effect | 0.154 | [ 0.069, 0.230] Total Effect | -0.178 | [-0.325, -0.027]

m16 <- brm(
  bf(ptgi_total_score ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m16)

summary(m16)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: ptgi_total_score ~ class 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    43.19      1.33    40.51    45.70 1.00     3225     2321
class2       -4.36      1.76    -7.72    -0.76 1.00     3158     2300

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    23.85      0.63    22.65    25.19 1.00     3179     2416
alpha     0.39      0.65    -0.73     1.68 1.00     2880     2317

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

bayes_R2(m16)

      Estimate   Est.Error         Q2.5      Q97.5
R2 0.009530306 0.006634892 0.0003035324 0.02528022

m17 <- brm(
  bf(ptgi_total_score ~ is_rescue_worker),
  family = skew_normal(),
  data = all_items,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

summary(m17)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: ptgi_total_score ~ is_rescue_worker 
   Data: all_items (Number of observations: 1068) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept             50.37      1.33    47.76    52.97 1.00     3038     2433
is_rescue_workerSi    -9.65      1.60   -12.78    -6.50 1.00     2988     2273

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    23.46      0.52    22.48    24.55 1.00     3236     2596
alpha     0.08      0.47    -0.76     0.98 1.00     3008     2787

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

m19 <- brm(
  bf(rate_of_activity_num ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m19)

Using 10 posterior draws for ppc type 'dens_overlay' by default.